Momentum-space finite-size corrections for Quantum-Monte-Carlo calculations 
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Extended solids are frequently simulated as finite systems with periodic boundary conditions, 
which due to the long-range nature of the Coulomb interaction may lead to slowly decaying finite- 
size errors. In the case of Quantum-Monte-Carlo simulations, which are based on real space, both 
real-space and momentum-space solutions to this problem exist. Here, we describe a hybrid method 
which using real-space data models the spherically averaged structure factor in momentum space. 
We show that (i) by integration our hybrid method exactly maps onto the real-space model periodic 
Coulomb-interaction (MPC) method and (ii) therefore our method combines the best of both worlds 
(real-space and momentum-space). One can use known momentum- resolved behavior to improve 
convergence where MPC fails (e.g., at surface-like systems). In contrast to pure momentum-space 
methods, our method only deals with a simple single-valued function and, hence, better lends itself 
to interpolation with exact small-momentum data as no directional information is needed. By virtue 
of integration, the resulting finite-size corrections can be written as an addition to MPC. 

PACS numbers: 71.10.Ca, 71.15.-m, 73.20.-r 



I. INTRODUCTION 

The issue of finite-size corrections in Quantum-Monte- 
Carlo (QMC) calculations has recently been attracting 
considerable attention,--'^"^ As QMC calculations of solids 
need to be carried out within a supercell, the Coulomb 
interaction is typically replaced by the so-called Ewald in- 
teraction that is compatible with the supercell geometry.^ 
This, in turn, is equivalent to dealing with an infinite 
system with a periodically repeated exchange-correlation 
(xc) hole. In other words, using the Ewald interaction 
includes the spurious effective interaction of an electron 
with its periodically repeated xc hole. One solution to 
this drawback - the real-space solution - is to use QMC 
data at length scales smaller than the supercell size where 
QMC is expected to be accurate and substitute the miss- 
ing terms implicitly or explicitly. E.g., the Model Peri- 
odic Coulomb-interaction (MPC) method^ deals with the 
periodically repeated xc hole by using the bare Coulomb 
term within the supercell and assuming that beyond the 
supercell an electron only feels the Hartree potential with 
no further correlations. This represents a good approxi- 
mation in a bulk solid, where the xc hole decays rapidly. 
In contrast, Chiesa ali^ traced the Ewald error back to 
an integration error and then added back the missing 
contribution, yielding a momentum-space solution. 

In a recent paper," the spherically-averaged structure 
factor Sk was introduced in order to analyze and reduce 
Coulomb finite-size effects. This {Sk) represents a natu- 
ral quantity to study finite size errors: On the one hand, 
it is a simple one-dimensional function which via inte- 
gration yields the interaction energy and, on the other 
hand, it naturally orders the QMC data according to 



length scale. Although QMC can, in principle, model 
correlations well at short length scales, it is bound to fail 
at larger scales due to the fact that finite simulation cells 
must be used. 

So far we have only touched on the issue of finite-size 
corrections for the interaction energy. The information 
to evaluate this is contained in the spherical average of 
the diagonal terms of the two-particle density matrix 
(the two particle density) and hence Sk, which contains 
the same information, suffices. However, other quanti- 
ties may be computed using QMC which also suffer from 
finite-size errors. E.g. Chiesa aL— also deal with correc- 
tions to the finite-size errors in the kinetic energy but 
since the kinetic energy needs information that goes be- 
yond the diagonal terms of the density matrix such an 
analysis is beyond the scope of this paper which is based 
solely on the information in Sk ■ Similarly, finite-size cor- 
rections to the momentum distribution'^ cannot readily 
be done using just Sk- 

The aim of this paper is to show that modelling Sk with 
the MPC ansatz yields a momentum-resolved MPC. This 
allows for an intuitive analysis of finite-size errors and for 
an improvement of MPC whenever MPC fails, such as in 
surface-like systems. First of all, in Section II we discuss 
Sk in the context of QMC and apply the MPC ansatz, 
thereby showing that via integration over k our method 
exactly maps onto MPC. In Section III, we first look at 
a non-interacting uniform electron gas in order to (i) see 
how MPC breaks down in this case and (ii) learn how 
to overcome the limitation of MPC by using our hybrid 
analysis. We also report Diffusion Monte Carlo (DMC) 
calculations of the structure factor of an interacting uni- 
form electron gas, which yields an intuitive understand- 
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ing of the MPC in interacting systems. This then leads 
us to an easy way to understand the existing difficulties 
that arise in the case of a slab geometry,i2iS and we pro- 
pose a simple method to address this difficult problem. 
We round off the paper with a summary and conclusions. 
We use atomic units throughout [fi — = m^. = 1). 



II. MODELLING THE EXCHANGE 
CORRELATION HOLE IN MOMENTUM SPACE 

A. The interaction energy 

The interaction energy [/™* of an arbitrary system of 
many interacting electrons can be expressed as follows 



(1) 



where Uh represents the Hartree energy (which in the 
case of an infinite jcllium model is exactly cancelled by 
the presence of the positive background) 
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(2) 



and Uxc represents the so-called xc interaction energy 
corresponding to the attractive interaction between each 
electron and its own xc hole: 



, , nj;c(r, r') 

dr' ^ / . (3) 



Uxc = ^ I drn{r) 



Here, n(r) is the electron density and nxc{'r,'r') repre- 
sents the xc-hole density of an electron at r. For brevity, 
we shall also define a reduced electron density Uredij, r'), 
which represents the electron density at r' seen by a given 
electron at r in the presence of (hence, reduced by) its 
xc hole: 



nred(r, r') = n(r') + rixdr, r') = 



n^(r, r') 
n(r) 



(4) 



where n^(r, r') is the so-called two-particle densityJl. Note 
that na;c(r, r') < 0. 
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FIG. 1. The reduced electron density nreci(i',r') seen by an 
electron at r, as a function of u = |r — r' |. In the absence of ex- 
change and correlation, this would equal the electron density 
n(r') (thin solid line) which in the case of a uniform electron 
gas would be constant. In the presence of exchange and cor- 
relation, the true reduced electron density rircdijc, r') is of the 
form represented by the thick line marked "TRUE", which 
for a supercell geometry (of linear dimension L) becomes the 
reduced Ewald density represented by the thick dotted line 
marked "EWALD" . The MPC interaction energy (7'"' is ob- 
tained by using a cutoff reduced electron density (thick solid 
line marked "MPC" ) that coincides with (i) the reduced Ewald 
electron density within the simulation cell and (ii) the electron 
density n{r') beyond the simulation cell. 



the sum of two terms: The Hartree energy Uh calculated 
with the Ewald interaction, and the beyond-Hartree xc- 
energy Uxc calculated with a cutoff Coulomb interaction 
using the minimum image convention, i.e., translating co- 
ordinates of electron pairs such that r — r' lies within the 
simulation cell. Hence, the MPC interaction energy J7™* 
is obtained by simply replacing the true reduced electron 
density nre(i(r,r') of Eq. (|4]) by the MPC reduced elec- 
tron density n^'l^'^ {r,r') of the form displayed in Fig. I 
by the thick solid line marked "MPC" . 



C. Spherically averaged structure factor Sk 



B. Model periodic Coulomb interaction (MPC) 

The only periodic solution of Poisson's equation for a 
periodic array of charges is the so-called Ewald interac- 
tion, which is of the Coulomb form 1 /r only in the limit 
of an infinitely large simulation cell. However, while the 
Hartree energy is given correctly by this Ewald inter- 
action, the part of the electron-electron energy coming 
from the interaction of electrons with their own xc hole 
yields a spurious contribution that is due to the inter- 
action of an electron with its periodically repeated xc 
hole. This drawback was solved in Ref. \M by replacing 
the Ewald interaction by a model periodic Coulomb in- 
teraction that yields an interaction energy consisting of 



Starting with the xc-hole density nxc{j^, r') at r' around 
an electron at r, one finds the following momentum- 
resolved form of the xc interaction energy of Eq. ([3]):^'^ 

Uxc^^J {Sk - 1) dk, (5) 

where Sk is the spherical average of the diagonal struc- 
ture factor S'k.k':'' 

Sk^l + ^J dvn{v) J dr'^^^n,,(r,r'), (6) 

and u = |r — r'|. 

Equation ([S]) represents a general expression for the xc 
interaction energy, which is formally exact not only for 
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homogeneous media but also for an arbitrary inhomoge- 
neous many-electron system. 

D. Monte-Carlo (MC) sampling of the structure 
factor 

When performing MC sampling on a function /(r, r') 
such as the Coulomb potential l/|r — r'|, what we are 
actually calculating is 

^/(ri,rj)\ = J dvn{v) J dr' n\r,r')f{r,r'), 

I MC 

(7) 

where the minimum image convention is used. In con- 
junction with Eq. (U), this yields the MC sampling of 
the spherically averaged structure factor of Eq. ([5]) : 

I MC 

here being the Hartree contribution: 

5f = l|^^drdun(r)n(r-u)^^, (9) 

where again the minimum image convention is being em- 
ployed, i.e., the integrals are carried out over the simula- 
tion cell (SC). The densities n(r) here are standard MC 
electron densities. 



E. MPC from Sk 

If one applies the MPC ansatz in the calculation of 
the spherically averaged structure factor Sk, by simply 
assuming (as in Fig. 1) that beyond the supercell corre- 
lations are only due to variations in the density, then 
Eq. ([5]) yields exactly the MPC interaction energy of 
Ref. [l|. This can be seen by introducing Eq. © into 
Eq. ^ and performing the k integration. Hence, the 
momentum-space based method of Ref. 4 reduces, un- 
der the MPC ansatz, to what we may call a momentum- 
resolved MPC. 

The MPC ansatz implies a finite extent of the xc hole, 
which in turn results in a quadratic behavior of S'fc as fc — >■ 
0. In the case of finite systems (e.g., atoms, molecules, 
and clusters) and bulk solids, this quadratic behavior of 
the structure factor is qualitatively correct because of 
the short range of the xc hole in those systems. Hence, in 
those systems, as long as the simulation cell is sufficiently 
large for the spherically averaged structure factor at the 
cutoff momentum kc ^ 1/L to already be in or close to 
the asymptotic low-fc behavior, the MPC ansatz yields 
accurate results: Since in those systems the true Sk is an 
essentially quadratic function as fc — ^ 0, constrained by 
Sk=o = 0, and the same constraints hold for 5*^^^*-^, only 
higher-order terms contribute to any residual error. 



There are some caveats, however. The non-interacting 
uniform electron gas as well as semi-infinite systems 
contain a linear contribution to 5*^ at fc — >■ (due to 
the presence of a xc hole that is not negligible even 
at large distances), which causes the failure of the 
MPC scheme. Nevertheless, in the framework of our 
momentum-resolved approach there is room for improve- 
ment over MPC, since one has flexibility to go beyond 
the MPC ansatz by replacing the low-fc structure factor 
Sk^^^ by its known correct value. 

III. THE MPC STRUCTURE FACTOR IN 
PRACTICE 

A. The non-interacting uniform electron gas 

The non-interacting (Hartree-Fock) uniform electron 
gas was dealt with extensively in Ref. i^. Here we briefly 
discuss this system, with the aim of now introducing a 
correction term. The exact Hartree-Fock (HF) structure 
factor Sk of a uniform electron gas, which is easily de- 
rived analytically^^ is shown in Fig. [5] (dashed-dotted 
red line) together with the result we obtain by using the 
MPC ansatz (solid black line). This figure clearly shows 
that the MPC result is in qualitative error at low fc, due 
to the fact that the actual HF structure factor does not 
exhibit a quadratic behavior in the limit as fc — >■ but 
a linear behavior instead. Our momentum-resolved tech- 
nique, however, has room for improvement, by going be- 
yond the MPC ansatz. 

Beyond a system-size dependent momentum cutoff kc 
(i.e., at fc > kc), the QMC electron correlation and hence 
structure factor are expected to be accurate. The mo- 
mentum cutoff fcc should be of the order of the inverse 
of the characteristic length L of the supercell. A sim- 
ple way to model a correction proceeds as follows: After 
choosing an appropriate cutoff fcc (vertical line of Fig. [5]), 
we model the corrected HF structure factor as a straight 
line between 5*^=0 = and and the MPC struc- 

ture factor as a quadratic curve between the same points. 
The correction to MPC is then given by the area between 
these two curves. The dotted blue line and the dashed 
green line of Fig. [5] illustrate this modelling. More real- 
istic models are of course feasible. E.g., one could take 
into account derivatives of 5*^ at fc = or fc = fcc, etc. 

Hence, as long as we are interested in the xc interac- 
tion energy Uxc, there is no need to actually evaluate the 
spherically averaged structure factor for all fc. In fact, if 
we assume that fcc is already in the asymptotic regime 
one can base the entire correction on the value of fcc and 
the asymptotic form of the structure factor. Let us as- 
sume the MPC structure factor below fcc is of the form 

S^''''' = Pk'+jk\ (10) 

while the true functional dependence (in the asymptotic 
region) ought to be 

Sk = ak , (II) 



MPC 
exact 




FIG. 2. (Color online) The spherically averaged structure fac- 
tor Sk of a non-interacting (Hartree-Fock) uniform electron 
gas with an electron-density parameter = 1 and 54 elec- 
trons in a face-centered-cubic (fee) simulation cell. The exact 
HF structure factor Sk is represented by a dashed-dotted red 
line and the MPC result by a solid black line. These two 
curves are also plotted in the inset, but now together with 
Eq. (1111) (nearly atop the exact HF structure factor), repre- 
sented by a blue dotted line, and the modelling of Eq. pop 
(nearly atop the MPC HF structure factor), represented by a 
green dashed line. The difference between the corresponding 
integrals represents the correction to the MPC HF data. The 
vertical line indicates the cutoff kc- 



a coming from the linear term in the HF Sk-^^ At the 
cutoff, both should of course coincide: 
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as should their derivatives: 



(12) 



(13) 



which yields /3 = 2a /k^ and 7 = —a/k^. 

In this model, the corrected structure factor Sk only 
differs from the MPC structure factor at momenta be- 
tween zero and k^, so the correction in Uxc can easily be 
evaluated using Eq. ([5]); it turns out to be: 



AUxc = u: 



corrected 



u. 



MPC 

xc 



— ak^ 
12 



(14) 



With kc ~ 1/L, we immediately get that the MPC error 
scales as ~ l/L^ ~ 1/^2/3 

We now also see how limiting the constraint of a finite 
cutoff when using an MPC like analysis actually is: The 
quadratic model HF structure factor of Eq. pUj) (sec the 
dashed green line of Fig. ^ and the actual MPC HF 
structure factor (solid black lines of Fig. ^ are nearly 
identical and quite different from the exact HF structure 
factor (dashed-dotted red line and dotted blue line of 
Fig. H]). This is despite the model MPC HF structure 
factor only using the asymptotic quadratic shape and 




FIG. 3. (Color online) DMC spherically averaged structure 
factor Sk of an interacting uniform electron gas of rs = 1 (in 
an fee simulation cell with A^=18, 54, 102, and 178, A'^ being 
the number of electrons), at the MPC level. The thick dotted 
blue line marked "asymptotic" represents the true structure 
factor of a uniform electron gas at — 0: Sk ~ /2u)p.— 
The vertical lines represent the cutoff kc — 1/L for 7V=18, 54, 
102, and 178. 



the value of the MPC HF structure factor at the cutoff 
where it ought to coincide with the exact HF structure 
factor. 

As pointed out above, so far this represents a crude way 
of devising the finite-size correction to the MPC result. 
More complex functional forms can be chosen, by using 
the structure factor or even its derivatives at fc = or fc 7^ 
to estimate the corresponding parameters. Analytic 
integration would then yield a more accurate correction 
term. 



B. The interacting uniform electron gas 

In Ref.lJ, Variational Monte Carlo (VMC) calculations 
of the structure factor Sk and the xc interaction energy 
Uxc of an interacting uniform electron gas were reported. 
Here we report DMC calculations of these quantities, as 
obtained by following the Hellmann-Feynman sampling 
introduced in Ref. [ij- 

Figure [3] exhibits the DMC structure factor Sk of an 
interacting uniform electron gas of Ts = 1, as obtained 
by following the MPC ansatz for various values of N (the 
number of electrons in the simulation cell). As pointed 
out in Ref. 4, the true interacting structure factor Sk is 
quadratic at A; ^ 0.^3 Fig. [3] shows that our MPC DMC 
structure factor nicely reproduces this \ow-k limit (thick 
dotted blue Hne marked "asymptotic"). The MPC DMC 
structure factor of Fig. [3] could be improved by using at 
low wavevectors the structure factor that one can obtain 
numerically in the random-phase approximation (RPA), 
which is known to be accurate in the low-fc regime. 
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C. Semi-infinite electron gas 

It is well known that in the case of a semi-infinite elec- 
tron gas the correct behavior of the structure factor Sk at 
low k consists of a linear term accounting for the surface 
contribution that is augmented by the usual quadratic 
bulk term. For a finite simulation cell of width and 
surface area A, one finds (in the long-wavelength limit 



where now 



and 



P = l/2wp. 



(15) 



(16) 



(17) 



Here, /3 represents the usual bulk term^i^ and a follows 
from Eqs. (2.13) and (3.34) of Ref.[il, = being 
the surface-plasmon energy. As the system gets larger 
{Lz — >■ c») the relative surface contribution shrinks, as 
expected. Nonetheless, surface energies, which have no 
contribution from the bulk part of the structure factor, 
are entirely dominated by the linear surface-contribution 
ak. 

The MFC ansatz will never yield the linear term of 
Eq. p^ . In order to solve this shortcoming, one can 
perform a simple analysis similar to the one leading 
to Eq. p4p . now employing the linear surface term of 
Eq. (|16p. This would yield an MFC error that scales as 
^ 1/V . However, for this to make sense one needs a well 
defined cutoff fcj, below which the quadratic MFC behav- 
ior should be replaced by the correct linear behavior of 
Eq. (fTS]) and above which the MFC structure factor is 
essentially exact. Nevertheless, in the case of realistic 
calculations^ of surface energies where the surface area 
needs to be varied for a fixed slab width and vice versa, 
finding a well-defined cutoff k^ might not be possible, so 
that one would need to explore more complex functional 
forms of the structure factor involving the MFC structure 
factor itself and possibly its derivatives at one or more 
values of k, such as A; = kc- In conjunction with the 
known surface contribution a of Eq. (|16|) . such schemes 
could correct significantly the existing MFC surface cal- 
culations. 



D. Implicit correction to the Coulomb kernel 

An interesting observation is that any correction to the 
structure factor implies a correction to the Coulomb 
potential V'^ = 1/u: For < /c < fcc , S^'^'-^ is given by 
evaluating the kernel sin(fczi)/fcu entering Eq. In- 
stead, sampling the quantity 



gcor. _ I ^^ sin(/fcu) 



S^'PC -I) ku 



(18) 



gives S^°^', the finite-size corrected structure factor. This 
will differ from S¥^'~^ only between k = Q and k — k^, 



the fraction in Eq. ()T8|) is therefore well-defined. But 
the integral of the kernel, (2/7r) " sm{ku)/ku = 1/u, is 
precisely the implied Coulomb potential V'^ . Changing 
the kernel thus changes the implicit Coulomb potential: 



dk 



Qcor. 



Sl'P^\ sin(fcii) 



gMPC _ 1 



ku 



(19) 



This is a short-range correction in the sense that when 
performing the MC run u is effectively restricted to 
within the simulation cell. On the other hand, the k 
in the integral is smaller than k^.^ which is inversely pro- 
portional to the simulation-cell size, and so ku entering 
the sin function in Eq. (fT9)) is generally below 1. 

In the case of the model of Eqs. PH)) and pT|) we get, 
for example. 



Ay" — 



dk 



2k 



-1 



fe3 



sm{ku). (20) 



Sk must then also be evaluated using AV'^\ as it is used 
in Sk via rixc- In contrast, Uh continues to be the stan- 
dard Ewald Hartree energy. 



IV. SUMMARY AND CONCLUSIONS 

First of all, we have demonstrated that modelling the 
spherically averaged structure factor Sk with the MFC 
ansatz (momentum-resolved MFC) yields exactly, after 
integration, the MFC interaction energy of Ref. |lj. This 
allows us to see explicitly that in the case of solids and 
finite systems MFC improves convergence considerably 
(over the more traditional Ewald scheme), due to the fact 
that for such systems MFC yields the correct quadratic 
behavior of S'fe as /c — ^ 0. We also see explicitly how 
the MFC ansatz breaks down in the case of the non- 
interacting uniform electron gas and, in general, in the 
case of all systems exhibiting a similar pathology, i.e., 
an extended xc hole (at surfaces, for example), where the 
leading term of 5^ at = is proportional to k. 

As the explicit k dependence of ^fc at low wave vectors 
can usually be known, one can look at QMC systems at 
different length scales separately enabling us to analyze 
the xc interaction energy at different length scales and to 
derive a correction term. On integration, this term yields 
a correction to QMC calculations that are based on the 
model periodic-Coulomb interaction MFC. 
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